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Dispersion properties of electromagnetic crystals formed by small uniaxial resonant scatterers 
(magnetic or electric) are studied using the local field approach. The goal of the study is to deter- 
mine the conditions under which the homogenization of such crystals can be made. Therefore the 
consideration is limited by the frequency region where the wavelength in the host medium is larger 
than the lattice periods. It is demonstrated that together with known restriction for the homoge- 
nization related with the large values of the material parameters there is an additional restriction 
related with their small absolute values. From the other hand, the homogenization becomes allowed 
in both cases of large and small material parameters for special directions of propagation. Two 
unusual effects inherent to the crystals under consideration are revealed: flat isofrequency contour 
which allows subwavelength imaging using canalization regime and birefringence of extraordinary 
modes which can be used for beam splitting. 

PACS numbers: 78.20.Ci, 42.70.Qs, 42.25.Gy 



INTRODUCTION AND PROBLEM 
FORMULATION 



The problem of homogenization of bulk arrays of small 
scatterers operating in the applied field as dipoles (el- 
electric or magnetic) has a long history. One can recall 
here the classical works of Lorentz, Madelung, Ewald and 
Oseen. However, in what concerns the homogenization 
of arrays of small resonant scatterers these classical re- 
sults were revised in 1970-s taking into account the pos- 
sible shortening the propagating wave at the resonance 
and the strong mutual coupling of resonant particles. It 
was done in the seminal work by Sipe and Kranendonk 
[1]. In 1990-s the interest to this problem was renewed 
by extensive studies of bianisotropic composites (see e.g. 
in [2] and [3]). The metal bianisotropic particles (chi- 
ral particles and omega particles) have small resonant 
size at microwaves due to their complex shape (they in- 
clude a wire ring and straight wire portions). However, 
the known studies of their homogenization are mainly re- 
ferred to the non-regular arrays. This can be explained 
by specific applications of microwave bianisotropic com- 
posites (as absorbers or antiradar coverings). The works 
like [4] concerning the regular bianisotropic lattices do 
not consider effects of particles resonance. Briefly, the 
homogenization problem for resonant scatterers has not 
been studied enough. However, it is becoming very im- 
portant now due to the following reasons. 

The first one is the rapid development of nano- 
technologies. It becomes possible to prepare the lattices 
of metallic nano-particles operating at the frequencies 
rather close to that of the plasmon resonance of the indi- 
vidual particle. Recently, a significant amount of works 
has been devoted to ID arrays (chains) of silver or gold 



particles which were found prospective for subwavelength 
guiding the light (see e.g. [5] and the list of references of 
this work). It is evident that the 2D and 3D lattices of 
metal nano-particles provide potentially even more broad 
scope of optical applications than the chains. If the ho- 
mogenization of a 3D lattice is possible then one can use 
the basic knowledge on the continuous media and apply 
it to the lattices. This approach can be rather instructive 
and we demonstrate below its example. In the present 
paper we study the case of microwave scatterers, but this 
is only an illustration of the theory. Similar results can 
be obtained for the optical range, too. The electric scat- 
terers of small resonant size are already known in optics, 
and the possibility to create the small resonant scatterers 
with magnetic properties was recently shown in work [6] . 

The second motivation of the present research is re- 
lated with the intensive studies of the so-called left- 
handed media [7] . The left-handed medium (LHM) is an 
effective continuous medium with simultaneously nega- 
tive permittivity and permeability. The all-angle nega- 
tive refraction and backward waves are inherent to such 
media. The interest to these artificial materials was 
evoked by seminal work of Pendry [8] indicating the op- 
portunity of the subwavelength imaging using a slab of 
LHM. The most loud realization of LHM is a uniaxial 
version of this medium studied in works [9], [10], [11] and 
others. This structure is composed from two compo- 
nents playing the roles of the building blocks. The first 
block (responsible for negative permittivity) is a wire 
medium [12-14] and the second block (responsible for 
negative permeability) is a lattice of the split-ring res- 
onators (SRR:s) [15]. The SRR:s are small magnetic scat- 
terers experiencing a two-time derivative Lorentz-type 
resonance. As a result, the permeability of the struc- 



ture can take negative values within the resonance band 
of SRR:s. This structure operates, however, as a LHM 
only in the plane orthogonal to wires (and for an only po- 
larization of the wave) . The reason of that is the strong 
spatial dispersion inherent to wire media at all frequen- 
cies [16]. This effect makes the axial component of the 
wire medium permittivity depending on the propagation 
direction. Only for the waves propagating in the orthog- 
onal plane this permittivity component is definitely neg- 
ative until the so-called plasma frequency and the struc- 
ture suggested in [9] can be treated as a LHM only in 
this special case. 

In order to obtain a variant of LHM operating in three 
dimensions some attempts to use small resonant electric 
scatterers together with magnetic ones [17] as well as the 
bianisotropic scatterers [18-20] were made. The samples 
of LHM obtained in [18, 20] demonstrate high losses in 
the LHM regime and this makes the known variants of 
isotropic LHM not very interesting. 

However, if the goal is to observe negative refraction 
and backward waves, and to obtain subwavelength im- 
ages in three dimensions, then the isotropic LHM is 
not an only solution. These effects can be obtained in 
anisotropic structures, too. And not only at high fre- 
quencies. The so-called indefinite media (in which the 
principal components of permittivity and permeability 
tensors have different signs) were studied in works [21- 
23]. These media offer variety of effects including neg- 
ative refraction, backward wave effect, near field focus- 
ing, high-impedance surface reflection, etc. Anisotropy 
of the media introduces additional freedom in manipula- 
tion by its dispersion and reflection properties [24] . Even 
a uniaxial media with negative permittivity along its axis 
allows to observe effects of negative refraction and back- 
ward wave with respect to the interface [25]. The the- 
oretical results [21-23] do not prove that the structure 
composed by wire medium and SRR:s will demonstrate 
these effects in practice. On the contrary, from [16] and 
[25] it is evident that these effects (which should exist in 
a continuous uniaxial medium with negative axial per- 
mittivity) are absent in wire media. In the same time, 
a lattice of uniaxial electric scatterers oriented in par- 
allel allows to obtain the negative axial permittivity for 
all directions of propagation (i.e. without spatial disper- 
sion) . If such a lattice substitutes the wire medium in the 
structure reported in [9, 10] then the effects predicted in 
[21-23] for continuous indefinite media can be obtained in 
practice. This is the second reason of the present study. 

In the current paper dispersion properties of electro- 
magnetic crystals formed by orthorhombic lattices of uni- 
axial magnetic or electric scatterers are studied. The 
orientation of scatterers along one of the crystal axes is 
considered. Geometry of the structure is presented in 
Fig. 1. As an example of magnetic scatterers we have 
chosen the SRR:s [9, 10, 15] (see Fig. 2). The electric 
dipoles are represented by the short inductively loaded 
wires (ILW) [17] (see Fig. 3). 

An analytical model based on dipole approximation 
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FIG. 1: Geometry of an electromagnetic crystal. The arrows 
show directions of dipole moments of the uniaxial scatterers. 



and local field approach is introduced. The dipole 
approximation (magnetic dipoles describing SRR:s and 
electric dipoles describing ILW:s) restricts the dimensions 
of inclusions to be much smaller than wavelength in the 
host media. The local field approach allows to take into 
account the dipolar interactions between scatterers ex- 
actly. It makes possible accurate studies of lattice reso- 
nances. The results allows to examine when the structure 
corresponds to its homogenized model of local uniaxial 
media and when not. 

It is well known that the lattices of resonant scatter- 
ers (though they do not exhibit the spatial dispersion at 
all frequencies unlike wire media) can experience spatial 
dispersion at low frequencies as compared to the spatial 
resonance of the lattice. This is the case when the wave- 
length in the medium becomes comparable with lattice 
period [1]. This results on resonance stop-band [1] and 
on appearance of complex modes within it [19, 26]. This 
makes the homogenization impossible within a sub-band 
belonging to the resonance band. In the present work 
we do not pay attention to the complex modes. The 
comparison of the original lattice and its homogenized 
model is made using the technique of isofrequency con- 
tours. Such an approach allows to check correspondence 
between properties of the structures under consideration 
and their homogenization models. 

Uniaxial media with negative permittivity (or perme- 
ability) along its axis and positive permittivity (or per- 
meability) in the transversal plane has the isofrequency 
contours of hyperbolic form [21-24]. These isofrequen- 
cies correspond to the negative refraction [22, 23]. If 
both original lattice and its homogenized model possess 
such isofrequencies then they both possess the negative 
refraction. More generally, if the homogenized model of 
the lattice keeps (at least approximately) the same isofre- 
quency contours, then the homogenization is allowed. If 
the homogenization dramatically change them the ho- 
mogenization is forbidden. This is the main idea of our 
approach. 
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II. MODELS OF INDIVIDUAL SCATTERERS 

The geometries of the SRRs and ILW are presented 
in Fig. 2 and Fig. 3, respectively. Since the dipoles 
moments of all scatterers are directed along x (see Fig. 
1) an individual scatterer can be characterized by scalar 
polarizability a relating the dipole moment with the local 
field (external field applied to a scatterer). 



A. Split-Ring Resonators 

The SRR considered in [9, 10, 15] is a pair of two copla- 
nar broken rings (see Fig. 2). Since the two loops are not 




FIG. 2: Geometry of a Split-Ring-Resonator 

identical the analytical models for such SRR:s are rather 
cumbersome [27, 28]. In fact, such SRR can not be de- 
scribed as a purely magnetic scatterer, because it ex- 
hibits bianisotropic properties and has resonant electric 
polarizability [27, 28] (see also discussion in [29]). How- 
ever, the electric polarizability and bianisotropy of SRR 
is out of the scope of this paper. We neglect these effects 
and consider an ordinary SRR as a magnetic scatterer. 
The analytical expressions for the magnetic polarizabil- 
ity a{ui) of SRRs with geometry plotted in Fig. 2 were 
derived and validated in [28]. The final result reads as 
follows: 



a(uj) — 
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where ujq is the resonant frequency of magnetic polariz- 
ability: 
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L is inductance of the ring (we assume that both rings 
have the same inductance): 



L = fi r 
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M is mutual inductance of the two rings: 



M = jji r 



(l-e)log(-) + C 
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C r is the effective capacitance of the SRR: 
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r is the radiation reaction factor: 
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r is the inner radius of the inner ring, w is the width of 
the rings, d is distance between the edges of the rings 
(see Fig. 2), e and no are permittivity and permeability 
of the host media, and k = io^/eofiQ is the wave number 
of the host medium. The presented formulae are valid 
within the frame of the following approximations: w,d <C 
r and the splits of the rings are large enough compared 
to d. Also, we assume that SRR is formed by ideally 
conducting rings (no dissipation losses). 

The magnetic polarizability (1) takes into account 
the radiation losses and satisfies to the basic Sipe- 
Kranendonk condition [1, 26, 30] which in the present 
case has the following form: 
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In the following analysis we operate with the inverse 
polarizability a _1 (w), thus, we rewrite (1) in the follow- 
ing form: 



a _1 (a;) = A- 



1 +j 



k 3 



(3) 



B. Inductively Loaded Short Wires 

A typical resonant electric scatterer is an inductively 
loaded short wire, as shown in Fig. 3. 
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FIG. 3: Geometry of the inductively loaded wire dipole 

The electric polarizability a e of loaded wires following 
the known model [17] has the form: 
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wire = n l £ o/ l°g(2V r o) i s the capacitance of the 



where C. 

wire, uiq = v 7 LC W i TC is the resonant frequency, L is the 
inductance of the load, I is the half length of the wire 
and ro is the wire radius. 

It is clear, that at the frequencies near the resonance 
the polarizability of LSW has the form 
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with A e = l 2 C w i YC , which is similar to (3). Moreover, if 
A e /£o = A/ no then using duality principle the magnetic 
dipole with polarizability a (3) can be transformed to the 
electric dipole with polarizability a e (3), and vice versa. 
This means that it is enough to consider only one type of 
resonant scatterers. In the present paper we have chosen 
magnetic ones to be principal. The electric scatterers 
can be easily obtained using duality principle from the 
magnetic scatterers with A = /ioA e /co- 



Ill. HOMOGENEOUS MEDIA APPROACH 

Let us consider an orthorhombic lattice with periods 
ax b x c formed by magnetic uniaxial scatterers directed 
along x (see Fig. 1) and described by polarizability (1). 
For electric scatterers (ILW:s) the problem is dual to the 
present one. In the long wavelength limit the lattices 
of scatterers are usually described as homogeneous me- 
dia with certain material parameters. The lattice under 
study can be modelled as a resonant uniaxial magnetic. 
The permeability of such a magnetic is a dyadic of the 
form: 

~p = ^x x + Mo (yoyo + z z ). 

The permeability [i (x-component of the tensor) can be 
calculated though the individual polarizability of a single 
scatterer using the Clausius-Mossotti formula [31]: 



M = Mo 
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(6) 



where V = abc is a volume of the lattice elementary 
cell and C s (a, b, c) is the static interaction constant of 
the lattice. The following expression for this interaction 
constant is available in [31], p. 758: 
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where Ko(x) is the modified Bessel function of 3d kind 
(the McDonald function) . In the case of a cubical lattice 
a = b = c the interaction constant is equal to the classical 
value C s = 1/(3V). 

Notice, that the radiation losses contribution in expres- 
sion (3) should be skipped while substituting into formula 
(6). This makes permeability purely real number as it 
should be for lossless regular arrays [1, 26]. This manip- 
ulation is based on the fact that the far-field radiation 
of the single scatterer is compensated by the electromag- 
netic interaction in a regular three-dimensional array, so 
that there is no radiation losses for the wave propagat- 
ing in the lattice. The mathematical proof of this fact 



for the general dimensions of lattice is presented in the 
Appendix. 

The typical dependence of magnetic permeability /i on 
frequency is presented in Fig. 4 for cubic lattice (a = b = 
c) of SRR:s with parameters chosen so that A = O.l^oa 3 
and luq = l/(a^/eo7*o) ■ The resonant frequency shift from 
ka = 1 to ka — 0.984 is clearly observed. While ka < 
0.984 the structure is paramagnetic (fx > 1). For ka 
within [0.984, 1.0352] range the permeability is negative 
(fi < 0). For ka > 1.0352 the medium is diamagnetic 
(0 <n< 1). 
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FIG. 4: Dependencies of relative permeability fx/^o and nor- 
malized polarizability a/(fioa [i ) vs. normalized frequency ka 
for cubic lattice (a — b — c) of SRR:s with A = 0.1/xoa 3 and 
loo = l/(a^/e /Lto). 

The dispersion equation for the uniaxial magnetic 
medium has the following form (see e.g. [24, 25, 31]): 



2 1 a 



n(k 2 



(8) 



Thus, the isofrequency surfaces for such material have 
form of a spheroid if /i > (the spheroid is prolate for /i < 
1 and oblate for fi > 1) or a hyperboloid if \i < 0. Both 
types of isofrequency surfaces have symmetry axis along 
OX. The media is isotropic in the YZ plane, and we can 
restrict our consideration by the XY plane without loss 
of generality. The typical isofrequency contours in this 
plane are shown in Figs. 5 and 6. The magnetic under 
consideration has the same parameters as in Fig. 4. The 
ranges of wave vector components q x and q y are restricted 
by ±7r/a and ±tt/6, respectively, having in mind that the 
exact dispersion diagram of the lattice corresponds to the 
first Brillouin zone and we will compare the homogenized 
model with the exact theory. 

While the frequency is below the resonance (ka < 
0.984) the isofrequency contour has the form of an el- 
lipse prolate along OY (/x > 1). For frequencies above 
the resonance but less than the frequency at which the 
permeability turns to zero (0.984 < ka < 1.0352) the 
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FIG. 5: Isofrequency contours in XY plane for uniaxial mag- 
netic with permeability as in Fig. 4 for frequencies near the 
resonance of the permeability. 
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FIG. 6: Isofrequency contours in XY plane for uniaxial mag- 
netic with permeability as in Fig. 4 for the frequencies where 
the permeability is close to zero. 



even at the frequency so that /i — > oo if one restricts by a 
special case of propagation. In general, it turns out that 
at all frequencies there are special cases of propagation 
for which the homogenization is allowed! 

The hyperbolic form of isofrequency contour is a 
unique feature inherent to resonant uniaxial magnetics 
[21-25]. It allows to achieve negative refraction at all in- 
cident angles for p-polarization in the case if the interface 
is normal to the optical axis. In the case of resonant uni- 
axial dielectric the same effect happens for s-polarization. 
If the uniaxial medium is two-component and has both 
negative axial permittivity and permeability, the nega- 
tive refraction should be observed for both p- and s- po- 
larizations [21-25]. 

Below we will compare Figs. 5 and 6 with those cal- 
culated for an original lattice of SRR:s using exact ap- 
proach. 



IV. DISPERSION EQUATION FOR 
ELECTROMAGNETIC CRYSTALS FORMED BY 
UNIAXIAL SCATTERERS 

Following the local field approach the dipole moment 
M of a zero numbered scatterer is determined by the 
magnetic field Hi oc . acting to this scatterer: M = aHf oc , 
where Hf oc = (Hi oc . ■ xrj). This local field is a sum of the 
magnetic fields H m nj ; produced at the coordinate origin 
by all other scatterers with indexes (to, n, I) ^ (0, 0, 0): 



H 



loc. 



(m,n,J)#(0,0,0) 



(9) 



The magnetic field produced by a single scatterer with 
index (to, n, I) is given by dyadic Green's function G(R): 



i = Mo G(R m ,„,;) M 



(10) 



where 



G(R) = ( 



k 2 I + VV ) 



-jkR 



) AttR 



isofrequency contours are hyperbolas (/i < 0), see Fig. 
4). If the frequency is above the frequency at which 
the permeability passes zero (ka > 1.0352) then the 
isofrequency contour becomes an ellipse oblate along OY 
(0 < /.t < 1). All the isofrequency contours pass through 
points q x = ±k. In particular, all the ellipses have the 
same semi-axes along x (equal to k). Notice, that the 
solution q x — k corresponds also to the arbitrary values 
of q y , q z if [i — > co. This implies the propagation of all 
waves along the optical axis x with same phase velocity 
which is equal to that of host medium. Strictly speaking, 
for an infinite value as well as for finite large values of /j, 
the homogenization is forbidden. But we will show using 
the local field method that the homogenization is allowed 



We consider uniaxial scatterers oriented along Xo di- 
rection, so it is enough to use only XoXo component of 
dyadic Green's function. So, we replace (10) by the scalar 
expression: 



where 
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To study eigenmodes of the system we introduce the 
phase distribution of dipole moments determined by the 
unknown wave vector q = (q x ,qy,qz) T as follows: 



M„ 



(12) 



Collecting together expressions (9), (11) and (12) we 
obtain dispersion equation relating the wave vector q 
with the frequency u>: 



M = afi 



(m,n, 0/(0, 0,0) 



G(R m , n ,i)ATe- 
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It can be rewritten in a more appropriate form: 

[ M o^ 1 M-C(A;,q)]M = 0, (13) 

where 

C(k,q,a,b,c)= Yl G(R m , nd )e-^ am+q y bn +^ cl l 

(m.n.l)=£(0,0.0) 

(14) 

We call C as the dynamic interaction constant of the 
lattice using the analogy with the classical interaction 
constant from the theory of artificial dielectrics and mag- 
netics [31]. 

Dispersion equation (13) has two different types of so- 
lutions. The first ones are ordinary waves with zero 
dipole moments (M — 0). They are plain waves prop- 
agating in the host media which have zero component 
of magnetic field along direction of dipoles. They do 
not interact with lattice (do not excite magneto-dipole 
moments). The waves of second type are extraordinary 
waves. They excite magneto-dipole moments (M ^ 0) 
strongly interacting with each other. The dispersion 
equation for extraordinary modes transforms from (13) 
to 



fioa (w) — C(k, q, a, b, c) = 0. 



(15) 



The solution of this dispersion equation allows to study 
dispersion diagrams for the crystal under consideration. 
The main problem is the calculation of the dynamic inter- 
action constant C given by (14). This question is closely 
related with such concepts as static interaction constant 
(7) and the triply-periodic dyadic Green's function. The 
static interaction constant can be obtained from (14) by 
letting k = q x = q y = q z = and choosing appropri- 
ate order of summation for obtained conditionally con- 
vergent series [32]. The plane- wise summation method 
[32, 33] or Poisson summation formulae based technique 
[31, 34] are usually applied for calculation of the static 
interaction constant. The triply-periodic dyadic Green's 
function represents the field produced by a phased lattice 
of point dipoles. If the zero-numbered term is added to 
series (14) (simultaneously one should move the observa- 
tion point in (10) from the node of the lattice to avoid 
singularity) formula (14) will give a co-polarized compo- 
nent of the dyadic Green's function. The triply-periodic 
dyadic Green's function is usually evaluated with help 



of classical Ewald's method [35-37]. However, the other 
methods of summation (with improved convergence rate) 
exist as well [38, 39]. All the listed above methods can 
be applied for evaluation of dynamic interaction constant 
(14). The Ewald's method require an appropriate choice 
of a splitting parameter [40] , which is a sophisticated ma- 
nipulation. Also, it does not show the energy balance in 
the lattice (as well as other methods mentioned above). 
Therefore, we have chosen an original method of summa- 
tion. Our approach combines the plane-wise summation 
[32, 33] and the Poisson summation technique with sin- 
gularity cancellation [31]. The details of the evaluation 
of C which includes the energy balance condition as an 
intermediate step are presented in Appendix. 



V. DISPERSION PROPERTIES OF THE 
CRYSTAL 

The dispersion equation (15) with interaction constant 
C(k, q) given by formula (A37) from Appendix is solved 
numerically. The parameters of the structure are the 
same as those of the homogenized structure: cubic lattice 
(a = b = c) of SRR:s, A = 0.1/ioa 3 andt^o = l/( a \/ £ o/- t o)- 
The dispersion diagram for the crystal is presented in 
Fig. 7. The points Y = (0,0, 0) T , X = (7r/a,0,0) T , Y 
= (0,7r/a,0) T , L = (0, n/a, n/af, K = (tt/o, 0, n/a) T 
and R = (7r/a, 7r/a, 7r/a) T of the first Brilloun zone are 
illustrated at the sketch in the left bottom corner of the 
plot. The dotted lines represent dispersion curves for 
ordinary modes of the crystal which coincide with light 
lines. The incomplete resonant stop band for extraordi- 
nary modes (similar to discussed in [26]) is observed in 
vicinity of the resonance of individual inclusions. 



light lines 




FIG. 7: Dispersion diagram for cubic lattice (a = b 
SRR:s with A = 0.1/xoa 3 and ujq = l/(a^/£oMo)- 



R 

c) of 



The dispersion curve for (010) direction (branch TY 
in Fig. 7) is shown in Fig. 8 for comparison with result 
predicted by homogenization model. For this direction of 
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FIG. 8: Dispersion curve for (010) direction. Exact solution 
(solid line), prediction by homogenization model (dashed line) 
and light line (dotted line). 
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FIG. 9: Isofrequency contours in YZ plane for frequencies 
near of the bottom edge of stop band. 



propagation the agreement with the homogenized model 
is fine except narrow frequency range ka ss 0.98. This re- 
gion is not visible in Fig. 8 but it is clear from Fig. 7 that 
this is the lower edge of the stop-band for waves propagat- 
ing in the transverse plane {YZ). This frequency range 
in the homogenization model corresponds to high propa- 
gation constant q y > n/a and high positive permeability 
H > it 2 / (ka) 2 . This means that the homogenization in 
the case /i > Tr 2 /(ka) 2 , strictly speaking, describes the 
dispersion of the lattice in a wrong manner. This is an 
expected result which corresponds to the known predic- 
tions of the classical theory [1]. Below we will consider 
this frequency range in details. 

The frequency band 0.9803 < ka < 1.044 corresponds 
to the negative axial permittivity of the homogenized 
model of the lattice. Negative axial permittivity means 
the imaginary propagation constant for the transverse 
plane and this result nicely corresponds to the stop-band 
for the YZ plane predicted by the exact theory. So, the 
homogenization within 0.9803 < ka < 1.044 is allowed. 

The isofrequency contours in YZ plane for frequen- 
cies near the bottom ka = 0.96 . . .0.9803 and top ka 
1 .()-", ... 1.10 edges of the transverse stop band are pre- 
sented in Fig. 9. The behavior of isofrequency contours 
shown in Fig. 9 is typical for general electromagnetic 
crystals at the frequencies near the stop band edges [41- 
43] . While the frequency is rather far below the stop band 
(ka < 0.979) the isofrequency contours have form of cir- 
cles and the agreement with the homogenized model is 
fine. The same behavior is observed above the stop band 
(ka > 1.044). The circles for ka > 1.044 are smaller 
than those for ka < 0.979 which nicely corresponds to 
the smaller effective permeability (see Fig. 4). However, 
within the narrow frequency range 0.979 < ka < 0.9803 
the isofrequency contours acquire a form which is differ- 



ent from a circle. This anisotropy in the transverse plane 
gives the evidence of spatial dispersion. Notice, that in 
this band in the lattice there are two evanescent modes 
whose wave vectors lie in the transverse plane (see also 
[26]). Strictly speaking, the crystal can not be homoge- 
nized at these frequencies. And these frequencies corre- 
spond to high positive /i of the homogenized lattice. It 
was already noticed above that it is the expected result. 
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FIG. 10: Isofrequency contours in XY plane for frequencies 
near the top edge of stop band. 

Significant disagreement between the exact solution 
and the result of homogenization was also obtained at 
the frequencies near the top edge of the stop band. The 
isofrequency contours in XY plane for this frequency 
range are presented in Fig. 10. They dramatically differs 
from prediction given by homogenization model shown 
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in Fig. 6. Following to homogenization approach, the 
isofrequency contours should have hyperbolic form at 
the frequencies corresponding to negative effective per- 
mittivity and elliptic one in the case of positive permit- 
tivity (see Fig. 6). The exact modeling reveals that 
this switching between hyperbolic and elliptic types of 
isofreqency contours happens in the different manner. 
While ka < 1.0435 the isofrequency contours has a form 
which is similar to hyperbolic one but they are already 
distorted. At the higher frequencies ka > 1.0435 the 
'hyperbolic' contours continue to distort, and simulta- 
neously the 'elliptic' contours (the second branch of the 
same isofrequency) appear in vicinity of T point. The 
'hyperbolic' contours pass through points q x = ±k while 
ka < 1.0455, but 'elliptic' ones do not. For ka > 1.046 
the situation changes to the opposite one. The 'ellip- 
tic' contours acquire fixed size along OX axis and starts 
to pass through points q x — ±fc. On the other hand, 
the 'hyperbolic' contours start to collapse around X-point 
and completely disappear for ka > 1.051. This way the 
'hyperbolic' contours transform to 'elliptic' ones passing 
through the regime where both types of contours co-exist 
at the same frequencies. At any frequency only one of 
these contours passes through points q x = ±k. 

Thus, in the region 1.043 < ka < 1.051 the homog- 
enization gives wrong results for the waves propagating 
in the XY plane because of the two-mode regime ob- 
served in original structure. This region corresponds to 
the small absolute values of fi (\fi\ < 0.2 in our case). 
Strictly speaking, the homogenization in the region of 
small turns out to be forbidden. In our opinion, this 
is a qualitatively new result. However, as it follows from 
Fig. 9 the homogenized model makes correct prediction 
for the waves propagating in the YZ plane in the band 
1.043 < ka < 1.051. One can conclude that the homoge- 
nization at these frequencies (forbidden in its strict mean- 
ing) is however allowed for a case of transversal propaga- 
tion. 

The described regime of the co-existence of 'hyperbolic' 
and 'elliptic' isofrequency contours at a fixed frequency 
means the bi-rcfringence for extraordinary modes and 
three-refringence in the case of the refraction (one or- 
dinary wave and two extraordinary ones). An extraordi- 
nary mode corresponding to the 'hyperbolic' contour re- 
fracts negatively, and the other one (corresponding to the 
'elliptic' contour) experiences positive refraction. This 
property can find different applications (beam splitting, 
etc). 



VI. CANALIZATION REGIME AND 
SUBWAVELENGTH IMAGING 

Above, we pointed out that near the bottom edge of the 
stop-band (frequencies corresponding to high positive /x) 
the homogenized model wrongly predicts the dispersion 
of the waves propagating in the YZ plane. Now let us 
show that the homogenized model gives the qualitatively 



correct predictions in this frequency region if considera- 
tion is restricted by the propagation in the XY plane. 
The isofrequency contours in XY plane for frequencies 
near bottom edge of stop band are presented in Fig. 11. 
The behavior of contours is in the good agreement with 
the predictions of the homogenized model (see Fig. 5). 
The difference is noticed only near the edges of the lowest 
Brillouin zone. So, the homogenization (forbidden in its 
strict meaning for ka ss 0.98) is still allowed for a case of 
oblique propagation with respect to the optical axis. 
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FIG. 11: Isofrequency contours in XY plane for frequencies 
near the bottom edge of stop band. 

At the frequencies near ka — 0.989 the isofrequency 
contours are practically flat. It means that all eigen- 
modes at such frequencies have the same axial component 
q x = ±k of the wave vector. Moreover, they all have the 
same group velocity (the group velocity is normal to the 
isofrequency contour). This makes the eigenmode to be 
the so-called transmission line mode like that of the wire 
medium [16]. For both lattice of uniaxial scatteres and 
wire medium this isofrequency corresponds to the infi- 
nite material parameter of the homogenized model. The 
difference is that in the present case the flat isofrequency 
contour exists at a single frequency ka — 0.989 in con- 
trast to wire medium which support transmission line 
modes in a very wide frequency range. 

The flat isofrequency contour we found can be used for 
the implementation of the so-called canalization regime 
described in our recent paper [43]. The similar regimes 
are called also as self-guiding [44], directed diffraction 
[45], self-collimation [46] and tunneling [47]. In [43] we 
have shown that not only all plane waves are collimated 
into a strictly parallel beam at the flat isofrequency. 
All evanescent waves impinging the medium at this fre- 
quency will be also transformed into the plane wave with 
q x = k transporting the energy along the optical axis. 
Therefore this regime allows to create subwavelength im- 
ages of the sources and transmit their near field to unre- 
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stricted distances. 

The canalization regime for a slab of the medium pos- 
sessing the flat isofrequency does not involve negative 
refraction and amplification of evanescent modes which 
are usually used for that purpose [8, 42, 48]. Its main 
feature is transformation of the spatial spectrum of the 
incident field into a collimated beam directed across the 
slab. All spatial harmonics of the source refract into such 
the eigenmodes at the front interface. These eigenmodcs 
all propagate normally to the interface with same veloc- 
ity and deliver the input distribution of electric field to 
the back interface. Their refraction at the back interface 
forms the image. The problem of the reflection from a 
slab (and inner reflections in the slab) can be solved using 
the Fabry-Perot resonance. The Fabry-Perot resonance 
holds for all incidence angles including the complex an- 
gles. The reason of that is simple: after the refraction all 
the incident waves acquire the same longitudinal compo- 
nent of the wave vector q x = k. Thus, in the canalization 
regime there is no image deterioration by the finite thick- 
ness of the lens (there are no waves travelling along the 
interfaces) . 

VII. CONCLUSION 

In the present paper we have studied dispersion prop- 
erties of the electromagnetic crystals formed by uniaxial 
resonant scatterers (magnetic and electric ones). The 
structures are modelled using the local field approach. 
The main tricky point of this theory is evaluating the dy- 
namic interaction constant of the lattice. This constant 
has been calculated using a special analytical method 
based on the plane-wise summation approach, Poisson 
summation formula, singularity cancellation technique 
and convergence acceleration of slowly convergent sc- 
ries. As a result, a transcendental dispersion equation 
has been obtained in the form suitable for rapid and ef- 
ficient numerical calculations. The comparison of exact 
solution provided by this equation with homogenization 



model allows to show that the structure, strictly speak- 
ing, can not be homogenized not only at the frequencies 
which correspond to the very high values of effective per- 
meability or permittivity (this was well known earlier) 
but also at the frequencies corresponding to small abso- 
lute values of them. 

However, if one is interested in a special cases of prop- 
agation then the homogenization can be allowed in both 
these frequency bands. For the propagation in the plane 
comprising the optical axis the homogenization is allowed 
in the region of large material parameters. For the prop- 
agation in the plane orthogonal to the optical axis the 
homogenization is allowed in the region of small material 
parameters. 

During this study we have found two interesting prop- 
erties of the crystals under consideration. At a single 
frequency near the bottom edge of stop band the isofre- 
quency contour is flat and this frequency corresponds 
to the infinite permeability or permittivity This fact 
makes possible to use the crystals for subwavelcngth 
imaging. The two-mode regime is observed at frequen- 
cies near the top edge of stop band. This corresponds to 
the bi-rcfringence for extraordinary waves and to three- 
rcfringencc of the incident wave in the general case of 
arbitrary polarization, which can be used for beam split- 
ting. 

The dispersion theory presented in this paper is a 
powerful tool for dispersion analysis of three-dimensional 
electromagnetic crystals. In the present form the theory 
is restricted to the case of simple (uniaxial) scatterers, 
but it can be extended to the case of electric or magnetic 
scatterers with arbitrary dyadic response. This will be 
done in our future publications. In this case it will be 
possible to develop an analytical theory for a lattice of 
the isotropic resonant scatterers (e.g. metallic spheres 
in the optical range) in more accurate manner than the 
known low-frequency approximations [31, 34, 35, 49] al- 
low to do. This can be actual for the optics of metal 
nano-particles which is developing fast. 
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APPENDIX A: EVALUATION OF DYNAMIC 
INTERACTION CONSTANT 



For calculation of the dynamic interaction constant 
C(k, q) (14) we apply a method based on plane-wise 
summation, Poisson summation formulae and singularity 
cancellation technique. This method was applied in [50] 
for calculation of the two-dimensional dynamic interac- 
tion constant for theory of doubly-periodic wire lattices. 
The series in (14) are divergent in classical meaning, but 
the physical reasoning of necessary type of summation 
is clear enough. Due to existence of losses in real space 
one should add infinitesimal imaginary part to the wave 
vector k of free space and tend it to zero in order to get 
the correct result. 

We split series (14) (remember that the zero term is 
excluded from summation) onto three parts: 
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FIG. 12: Splitting areas 

These parts are denoted as £-1,2,3 respectively, and 
C = C\ + C2 + C3. The splitting areas are shown in 
Fig. 12. The term C\ describes contribution into the lo- 
cal field from all plane grids which are parallel to the 
OXY plane except the grid located at this plane. The 
term C2 corresponds to the contribution of the dipolc lin- 
ear chains parallel to the OX axis and located at OXY 
plane except the chain located at this axis. The term 
C3 is the contribution from all dipoles of the chain lo- 
cated at the OX axis except the dipole located at the 
origin of coordinate system. Notice, that C2 + C3 gives 
the interaction constant of the planar grid. 

For evaluation of the term C\ it is possible to use the 
Poisson summation formula for double series which leads 
to the expression with rapidly (exponentially) convergent 
series. The term C2 can be calculated using the ordinary 
Poisson summation formula together with the singular- 
ity cancellation technique [31]. It is impossible to apply 
Poisson summation formula for evaluation of the term 
C3 since it contain non-complete series. Convergence of 
these series can be accelerated using the dominant part 
extraction [51]. 
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1. Contribution of parallel planar grids 

The double Fourier transformation of any function 
f(x,y) is defined as follows: 

+00 +00 

F(p,q)=-L x , v {f(x,v)}= J J f(x,y)e-^ x+ ^dxdy. 



The Poisson summation formula for double series has 
the following form: [31]: 
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The double Fourier transform of the Hertz potential of 
a dipole reads as: 
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where the sign of square root should be chosen so that 



Im (yp 2 + q 2 - k 2 j > 0. 

Applying the shift and differential properties of Fourier 
transformation to the Fourier-image of the Hertz poten- 
tial (A2) we obtain the transformation rule: 
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Using the Poisson summation formula (Al) together 
with (A3) we obtain the term C\ in the form: 
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where 



2im 



k 2 , 



fc(™") = -jy(^) -k 2 . 

Here we choose Im(\A) > 0, so that Im (k z mn ^ < 0. 
The representation (A4) can be treated as an expansion 



of the fields produced by parallel dipole grids in terms of 
the Floquet modes. The wave vectors of these modes are 

k (mn) = (fcM ; fc(") ; fc(™«)) T . 

The series by I index in (A4) are geometrical progres- 
sions and their summation can be made directly: 
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It allows to rewrite expression (A4) for the term C\ as 
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These series possess exponential convergence. It is clearly 
seen if the second factor of the term under sign of the sum 
in (A5) is rewritten as 
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This makes (A5) suitable for rapid numerical calcula- 
tions. 



2. Contribution of parallel chains from OXY-plane 

The ordinary Fourrier transformation has the form 

+00 

F(p)=L x {f(x)}= J f{x)e-™*dx. 

— CO 

Poisson's summation formula for single series reads: 
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The Fourrier transform for the Hertz's potential of a 
dipole is following: 
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(A7) 

Thus, applying shift and differential properties of Fur- 
rier transformation to the image of Hertz's potential (A7) 
we get the following transformation rule: 
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Using Poisson's summation formulae (A6) together 
with (A8) we obtain the term C 2 in the form: 



+00 



C 2 = -]T J2 ^f-K ( Pm \bn\)e-^. (A9) 
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If arguments of McDonald's functions in (A9) have 
nonzero real part then the series by index n have very 
good convergence, but if these are getting imaginary 
then McDonald's functions transform to Hankcl's func- 
tions and the mentioned series become slowly convergent. 
Therefore we separate the part of (A9) which has good 
convergence: 
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which has slow convergence should be calculated with 
help of special method. Note, that there is only finite 
number of indexes m such that Re(p m ) = 0. It means 
that in (All) the summation by index m includes only 
finite number of terms. For example, at the low frequency 
limit, when period a is large compared with wavelength 
in the host media, the equation Ke(p m ) = has only one 
solution m = if q x < k. 

We will calculate the sum of the scries (All) as the 
limit with z tending to zero: 
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Introducing the auxiliary parameter z makes possible to 
complement series (A12) by zeroth terms and then to use 
the Poisson summation formula by index n (see (A3) for 
necessary Fourier transform). The result is as follows: 
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The term K (p m \z\) in (A13) plays role of the zero 
term which is subtracted from the complete series (al- 
ready transformed using Poisson summation formula) in 
order to get series (A12) without zero term. This term 
contains singularity if z tends to zero. This singularity 
disappears in (A13) during substraction from the com- 
plete series which experience the same singularity. In or- 
der to cancel out these singularities analytically we apply 
the method of a dominant series. Namely, we split the 
series from (A13) onto dominant and correction parts: 
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The dominant series can be evaluated using the tabu- 
lated formula (see [31], Appendix): 
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The whole singularity is included into the dominant 
series. The correction series have no singularity when z 
tends to zero. Using this fact the formula (A13) can be 
rewritten as: 
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The logarithmic singularity occurring here is compen- 
sated by that arising from the term with the McDonald 
function. The small-argument expression for the McDon- 
ald function reads: 

K (a) -> - [7 + Iog(a/2)] , 

where 7 = 0.577 is Euler's constant. 

Thus, the value of the limit in (A15) is as follows: 
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The series in (A15) with index n G [—00, +00] except 
n = have convergence 1/n 2 . Such convergence rate 
makes calculations not rapid enough, the convergence 
can be accelerated by extraction of the dominant series. 
In order to get the convergence 1/n 4 it is enough to ex- 
tract series of order 1/n 2 and 1/n 3 : 
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where l m = 2q 2 — p m and we have taken into account, 
that 



n=l 
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Collecting the terms corresponding to +n and — n to- 
gether in (A15) we obtain: 
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The property /ci" 1 ' n \q x ,q y ) = ki m ' n \q x , -q y ) makes 
function 
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even with respect to variable q y b/(2Trn). It means that 
being expanded into Taylor series it will contain only even 
power terms. Thus, the transformed series (A18) being 
expanded as the series of the order 1/n will contain only 
odd-power terms. In (A18) we have already extracted the 
dominant series of the order 1/n 3 . So, we conclude that 
series (A18) have convergence of the order 1/n 5 which is 
better than it was estimated when we started to extract 
dominant series in (A17) 

Collecting the parts of the term Ci given by (A10), 
(A15), (A16) and (A18) we obtain the final formula for 
Ci possessing convergence which is appropriate for rapid 
and effective numerical calculations: 
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The first series in the right-hand side of (A21) (that 
containing the expression in prances) can be simplified 
up to 
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These series have convergence 1/m 4 which is convenient 
for rapid calculations. The other series in the right-hand 
side of (A21) can be evaluated in the closed form using 
the formula (A14): 
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3. Contribution of the line located at OX-axis 



The term C3 has form of the series [5, 52]: 
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These series have convergence which is not enough for 
effective direct numerical calculations. We will use con- 
vergence acceleration technique presented in [51] in order 
to evaluate these series. The dominant series can be ex- 
tracted in the following way: 
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After these manipulations the formula (A21) trans- 
forms as follows: 
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where 



t+ = 1 - e j ( k+q * )a , t- = 1 - eJ( fc -9-) a . 

The expression (A22) looks more cumbersome as com- 
pared to the initial formula (A20), but it is much more 
convenient for rapid calculations. The estimations show 
that in order to get accuracy of 0.01% one needs to take 
more than 200 terms in expression (A20) and only 10 
terms in (A22). 
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4. Energy conservation 



The imaginary part of formulae (A5) reads: 



In this subsection we evaluate the imaginary part of C 
and consider the problem of the energy balance in a ID 
array (chain) of dipoles, in a 2D array (grid) and in a 3D 
array (lattice). 

Let us return to formula (A20) and find the imaginary 
part of the interaction constant of the dipole chain: 
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To calculate these series we used the auxiliary formulas 
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(A25) 



These formulas can be easily obtained from the relation 
(A14) rewritten for the case a = js 
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where s' = 2ir{s/(2ir)} and we use notation {x} for frac- 
tional part of variable x. To derive (A24) and (A25) one 
should integrate (A26) by parameter s. 
Note, that real part of C3 contain series 
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which can be expressed in terms on second and third 
order repeating integrals of tangent function. These in- 
tegrals can not be evaluated in elementary functions, but 
they are suitable for numerical calculations. We prefer 
to use the acceleration technique leading to the result 
(A22) for evaluation of C3, rather than to use numerical 
integration. Some other recommendations for calculation 
of series (A27) can be found in [31] together with their 
expansions into Taylor series. 

After substitution of (A24) and (A25) into (A23) and 
some algebra the following compact form for Im(Cs) can 
be obtained: 
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It is easy to obtain the imaginary parts of C 2 and C\ . 
The imaginary part of formula (A19) reads as: 
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Collecting together (A28), (A29) and (A30) we obtain 
that 



Im(C) = 
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(A31) 



This relation makes dispersion equation (13) real valued 
for the case of propagating modes. 

Now, let us discuss the energy balance in the chain 
using the result (A28). If the dipoles are arranged in a 
periodical linear array x = am phased by wave vector 
with the x— component q x (as in [52]) then the structure 
radiates cylindrical waves. The number of these waves 
depends on the relation between the wavelength, chain 
period and phase constant q x . In the regime of the guided 
mode q x > k this number is zero since \k x m) \ > k for 
all m. Using the Sipe-Kranendonk condition (2) for the 
imaginary part of the polarizability's inverse value one 
can obtain a purely real valued dispersion equation for 
the guided mode in the chain [5, 52]: 

Aioa _1 (w) - C 3 (u), q x , a) = 

However, the arrangement of the dipoles into an array 
changes the radiation losses of the individual scatterers. 
The effective polarizability of the scatterer in the linear 
array becomes as follows: 

011 = (o^ 1 - ^C^) 1 . 

The Sipe-Kranendonk condition in the general case of 
radiated waves should be replaced by 
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The expression (A32) follows from the formulae (2) and 
(A28). This relation expresses the balance between the 
radiation losses of the individual scatterer of the chain 
and the contribution of the chain unit cell into the radi- 
ated waves (|fci m ' l | < k). 

Now, consider a 2D grid of dipoles located at the nodes 
with coordinates x = am and y = bn and phased by real 
q x and q y , respectively. The effective polarizability of a 
scatterer in this planar grid is 

a 2 = K 1 - Mo 1 C 2 y 1 = (a- 1 - ^{C 2 + C 3 )) _1 . 

(A33) 

The formulae (A33) and (A29) allow to formulate an 
analogue of the Sipe-Kranendonk condition for the planar 
grid: 
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The terms — p^/(4a/uo) corresponding to the cylindrical 
waves in (A32) are cancelled out by the respective terms 

from (A29) and replaced by terms p^/ (2ab/j,oki mn ^). The 
last ones correspond to the radiated plane-waves (Flo- 
quet harmonics with indexes (m, n) produced by the 

grid). The condition Im(fci" m ' ) ) = for the finite sum 
in (A34) is the radiation condition for these Floquet har- 
monics. Formula (A34) expresses the balance between 
the radiation losses of the dipole and the contribution of 
the grid unit cell into radiation. 

In the surface wave regime, when Im(fci m ™') ^ for all 
m,n, using the Sipe-Kranendonk condition (2) one can 
obtain a real valued dispersion equation for the surface 
wave propagating along the grid: 



fi a 1 (oj) - C2{co,q x ,q y ,a,b) = 0, 



where 



C2{uj,q x ,q y ,a,b) = C 2 {uj,q x ,q y ,a,b) + C 3 (ui,q x ,a). 

Finally, let us consider a 3D lattice with nodes x = 
am, y — bn and z — cl phased by real q x , q y and q z 
respectively. The effective polarizability of the scatterer 
in this lattice is 



a 3 = K 1 - Ci) 1 - (a^-Cy 



(A35) 



From (A30),(A35) and (A34) we easily obtain 

Im(a3 1 )=0. (A36) 

The terms p 2 n /(2abjj,ok z mn ^) in (A34) are cancelled by 
respective terms of (A30). Physically, it means that ra- 
diation losses of the scatterer in this lattice are zero. The 
lattice does not radiate power because it fills the whole 
space and the radiation losses of the single scatterer are 
compensated by the electromagnetic interaction in the 
lattice (as well as in the waveguide regimes of the chain 
and of the grid). 



7T7J 87T 3 n 3 



l m b 3 b / Hp 
1.202^V + - log 1 



87T 3 7T 



47ra 3 



(2jka + 3)m + 2 
2-~> m 3 (m + l)(m + 2) 

m=l x ' K ' 



4n + V +J 2 



e -jkam cog ( feam ) 



-(jka + 1) (4 logt+ + t 2 _ log* - + 2eJ ka cos(q x a)) 



-2jka (t+ \ogt+ + t- \ogt~) + (7 jka + 3) 



where we use following notations (introduced above and 
collected here): 
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The calculations using (A37) can be restricted to the 
real part only, because its imaginary part is predefined 
by (A31). The series in (A37) have excellent convergence 
that ensure very rapid numerical calculations. 



5. Final formula 



6. Low frequency limit case 



Collecting together results (A5), (A19), (A20) we ob- 
tain the final expression for the dynamic interaction con- 
stant: 
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It is useful to consider the low frequency limit (when 
k, q Xl q y and q z are small as compared with 1/a, 1/b and 
1 /c) and show that the result for C transits to the known 
one for this case. Following to definition (A20) for term 
C3 we conclude, that 
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The expression (A19) for C 2 reduces to 
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Note, that both C3 and C2 turn out to be independent 
on k and q. The formula (A5) for C\ splits into the two 
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terms: the first one which depend on k and q (where we 
have expanded trigonometric functions into Taylor series) 
and some additional constant: 
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+ C s {a,b,c), (A38) 



where C s (a 1 b,c) is static interaction constant (7), and 
obtain alternative representation for C s (a, b, c): 
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(A39) 

The static interaction constant expressed as (A39) is 
equivalent to (7). The expression (A39) can be obtained 
from (7) applying Poisson summation formula by index 
n and then by direct summation by index I (in the same 
manner as it was done above during evaluation of term 
Ci). Both formulae (A39) and (7) are extremely effective 
for rapid numerical calculations due to excellent conver- 
gence of the series. The difference between (7) and (A39) 
is that (7) contains triple series in contrast to (A39) which 
comprises only double ones. Noteworthy, that conver- 
gence of series in (7) is higher than in (A39). 

Formula (A38) being substituted into (13) reduces the 
dispersion equation for an electromagnetic crystal to the 
known dispersion equation of a continuous uniaxial mag- 
netic (8) with magnetic permittivity of the form (6). This 
fact is an important verification of introduced dispersion 
theory. 



